Covariates with Missing Values

Estimate Feb SWE using AR(1) with spatial correlation

Imagine that for our study we need an estimate of SWE for all sites. We will use the information from the sites with full data to estimate the missing SWE for other sites. We will use a MARSS model to use all the available data.

\[\begin{equation} \begin{bmatrix} x_1 \\ x_2 \\ \dots \\ x_{15} \end{bmatrix}_t = \begin{bmatrix} b&0&\dots&0 \\ 0&b&\dots&0 \\ \dots&\dots&\dots&\dots \\ 0&0&\dots&b \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \dots \\ x_{15} \end{bmatrix}_{t-1} + \begin{bmatrix} w_1 \\ w_2 \\ \dots \\ w_{15} \end{bmatrix}_{t} \\ \begin{bmatrix} y_1 \\ y_2 \\ \dots \\ y_{15} \end{bmatrix}_t = \begin{bmatrix} x_1 \\ x_2 \\ \dots \\ x_{15} \end{bmatrix}_t + \begin{bmatrix} a_1 \\ a_2 \\ \dots \\ a_{15} \end{bmatrix}_{t} + \begin{bmatrix} v_1 \\ v_2 \\ \dots \\ v_{15} \end{bmatrix}_t (\#eq:mssmiss-ar1) \end{equation}\]

We will use an unconstrained variance-covariance structure for \(\mathbf{w}\) and assume that \(\mathbf{v}\) is identical and independent and very low (SNOTEL instrument variability). The \(a_i\) determine the level of the \(x_i\).

We need our data to be in rows. We will use reshape2::acast().

dat.feb <- reshape2::acast(swe.feb, Station ~ Year, value.var = "SWE")

We set up the model for MARSS so that it is the same as @ref(eq:mssmiss-ar1). We will fix the measurement error to be small; we could use 0 but the fitting is more stable if we use a small variance instead. When estimating \(\mathbf{B}\), setting the initial value to be at \(t=1\) instead of \(t=0\) works better.

ns <- length(unique(swe.feb$Station))
B <- "diagonal and equal"
Q <- "unconstrained"
R <- diag(0.01, ns)
U <- "zero"
A <- "unequal"
x0 <- "unequal"
mod.list.ar1 = list(B = B, Q = Q, R = R, U = U, x0 = x0, A = A, 
    tinitx = 1)

Now we can fit a MARSS model and get estimates of the missing SWEs. Convergence is slow. We set \(\mathbf{a}\) equal to the mean of the time series to speed convergence.

library(MARSS)
m <- apply(dat.feb, 1, mean, na.rm = TRUE)
fit.ar1 <- MARSS(dat.feb, model = mod.list.ar1, control = list(maxit = 5000), 
    inits = list(A = matrix(m, ns, 1)))

The \(b\) estimate is ````.

Let’s plot the estimated SWEs for the missing years (Figure @ref(fig:mssmiss-snotelplotfits-ar1)). These estimates use all the information about the correlation with other sites and uses information about correlation with the prior and subsequent years. We will use the tidy() function from the broom package to get the estimated 95% confidence intervals for the estimated states. Notice that for some sites, CIs are low in early years as these sites are highly correlated with site for which there are data. In other sites, the uncertainty is high in early years because the sites with data in those years are not highly correlated.

(ref:mssmiss-snotelplotfits-ar1) Estimated SWEs from the expected value of the states \(\hat{x}\) conditioned on all the data. Note we could also use the expected value of the \(y\) conditioned on all the data.

fit <- fit.ar1
d <- augment(fit, interval = "confidence")
d$Year <- d$t + 1980
d$Station <- d$.rownames
p <- ggplot(data = d) + geom_line(aes(Year, .fitted)) + geom_ribbon(aes(x = Year, 
    ymin = .conf.low, ymax = .conf.up), linetype = 2, alpha = 0.5)
p <- p + geom_point(data = swe.feb, mapping = aes(x = Year, y = SWE))
p + facet_wrap(~Station) + xlab("") + ylab("SWE (demeaned)")

(ref:mssmiss-snotelplotfits-ar1)

If we were using these SWE as covariates in a site specific model, we could then use the estimates as our covariates, however this would not incorporate uncertainty. Alternatively we could use Equation @ref(eq:mssmiss-marsscovarx) and set the parameters for the covariate process to those estimated for our covariate-only model. This approach will incorporate the uncertainty in the SWE estimates in the early years for the sites with no data.

Note, we should do some cross-validation (fitting with data left out) to ensure that the estimated SWEs are well-matched to actual measurements. It would probably be best to do ‘leave-three’ out instead of ‘leave-one’ out since the estimates for time \(t\) uses information from \(t-1\) and \(t+1\) (if present).

Diagnostics

The state residuals have a tendency for negative autocorrelation at lag-1 (Figure @ref(fig:mssmiss-stateresids-ar1)).

(ref:mssmiss-stateresids-ar1) State residuals for the AR(1) model. Many stations for autocorrelation at lag-1.

fit <- fit.ar1
par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
apply(residuals(fit)$state.residuals[, 1:30], 1, acf)

(ref:mssmiss-stateresids-ar1)

Estimate Feb SWE using only correlation

Another approach is to treat the February data as temporally uncorrelated. The two longest time series (Paradise and Olallie Meadows) show minimal autocorrelation so we might decide to just use the correlation across stations for our estimates. In this case, the state of the missing SWE values at time \(t\) is the expected value conditioned on all the stations with data at time \(t\) given the estimated variance-covariance matrix \(\mathbf{Q}\).

\[\begin{equation} \begin{bmatrix} x_1 \\ x_2 \\ \dots \\ x_{15} \end{bmatrix}_t = \begin{bmatrix} w_1 \\ w_2 \\ \dots \\ w_{15} \end{bmatrix}_{t}, \, \begin{bmatrix} w_1 \\ w_2 \\ \dots \\ w_{15} \end{bmatrix}_{t} \sim \begin{bmatrix} \sigma_1&\zeta_{1,2}&\dots&\zeta_{1,15} \\ \zeta_{2,1}&\sigma_2&\dots&\zeta_{2,15} \\ \dots&\dots&\dots&\dots \\ \zeta_{15,1}&\zeta_{15,2}&\dots&\sigma_{15} \end{bmatrix} \\ \begin{bmatrix} y_1 \\ y_2 \\ \dots \\ y_{15} \end{bmatrix}_t = \begin{bmatrix} x_1 \\ x_2 \\ \dots \\ x_{15} \end{bmatrix}_t + \begin{bmatrix} a_1 \\ a_2 \\ \dots \\ a_{15} \end{bmatrix}_{t} + \begin{bmatrix} v_1 \\ v_2 \\ \dots \\ v_{15} \end{bmatrix}_t (\#eq:mssmiss-tempind) \end{equation}\] Again \(\mathbf{a}\) is the mean level in the time series. Note that the expected value of \(\mathbf{x}\) is zero if there are no data, so \(E(\mathbf{x}_0)=0\).

ns <- length(unique(swe.feb$Station))
B <- "zero"
Q <- "unconstrained"
R <- diag(0.01, ns)
U <- "zero"
A <- "unequal"
x0 <- "zero"
mod.list.corr = list(B = B, Q = Q, R = R, U = U, x0 = x0, A = A, 
    tinitx = 0)

Now we can fit a MARSS model and get estimates of the missing SWEs. Convergence is slow. We set \(\mathbf{a}\) equal to the mean of the time series to speed convergence.

m <- apply(dat.feb, 1, mean, na.rm = TRUE)
fit.corr <- MARSS(dat.feb, model = mod.list.corr, control = list(maxit = 5000), 
    inits = list(A = matrix(m, ns, 1)))

The estimated SWEs for the missing years uses the information about the correlation with other sites only.

fit <- fit.corr
d <- broom::augment(fit, interval = "confidence")
d$Year <- d$t + 1980
d$Station <- d$.rownames
p <- ggplot(data = d) + geom_line(aes(Year, .fitted)) + geom_ribbon(aes(x = Year, 
    ymin = .conf.low, ymax = .conf.up), linetype = 2, alpha = 0.5)
p <- p + geom_point(data = swe.feb, mapping = aes(x = Year, y = SWE))
p + facet_wrap(~Station) + xlab("") + ylab("SWE (demeaned)")

fit <- fit.corr
par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
apply(residuals(fit)$model.residuals[, 1:30], 1, acf)

Estimate Feb SWE using DFA

Another approach we might take is to model SWE using Dynamic Factor Analysis. Our model might take the following form with two factors, modeled as AR(1) processes. \(\mathbf{a}\) is the mean level of the time series.

\[ \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}_t = \begin{bmatrix} b_1&0\\0&b_2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}_{t-1} + \begin{bmatrix} w_1 \\ w_2 \end{bmatrix}_{t} \\ \begin{bmatrix} y_1 \\ y_2 \\ \dots \\ y_{15} \end{bmatrix}_t = \begin{bmatrix} z_{1,1}&0\\z_{2,1}&z_{2,2}\\ \dots\\z_{3,1}&z_{3,2} \end{bmatrix}\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}_t + \begin{bmatrix} a_1 \\ a_2 \\ \dots \\ a_{15} \end{bmatrix} + \begin{bmatrix} v_1 \\ v_2 \\ \dots \\ v_{15} \end{bmatrix}_t \]

The model is set up as follows:

ns <- dim(dat.feb)[1]
B <- matrix(list(0), 2, 2)
B[1, 1] <- "b1"
B[2, 2] <- "b2"
Q <- diag(1, 2)
R <- "diagonal and unequal"
U <- "zero"
x0 <- "zero"
Z <- matrix(list(0), ns, 2)
Z[1:(ns * 2)] <- c(paste0("z1", 1:ns), paste0("z2", 1:ns))
Z[1, 2] <- 0
A <- "unequal"
mod.list.dfa = list(B = B, Z = Z, Q = Q, R = R, U = U, A = A, 
    x0 = x0)

Now we can fit a MARSS model and get estimates of the missing SWEs. We pass in the initial value for \(\mathbf{a}\) as the mean level so it fits easier.

library(MARSS)
m <- apply(dat.feb, 1, mean, na.rm = TRUE)
fit.dfa <- MARSS(dat.feb, model = mod.list.dfa, control = list(maxit = 1000), 
    inits = list(A = matrix(m, ns, 1)))

Diagnostics

The state residuals are uncorrelated.

fit <- fit.dfa
par(mfrow = c(1, 2), mar = c(2, 2, 1, 1))
apply(residuals(fit)$state.residuals[, 1:30, drop = FALSE], 1, 
    acf)

As are the model residuals:

par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
apply(residuals(fit)$model.residual, 1, function(x) {
    acf(na.omit(x))
})

Modeling season across sites

We will create a state for the seasonal cycle and each station will have a scaled effect of that seasonal cycle. The observations will have the seasonal effect plus a mean and residuals (observation - season - mean) will be allowed to correlate across stations.

ns <- dim(dat.yr)[1]
B <- "zero"
Q <- matrix(1)
R <- "unconstrained"
U <- "zero"
x0 <- "zero"
Z <- matrix(paste0("z", 1:ns), ns, 1)
A <- "unequal"
mod.list.dfa = list(B = B, Z = Z, Q = Q, R = R, U = U, A = A, 
    x0 = x0)
C <- matrix(c("c1", "c2"), 1, 2)
c <- c.seas
mod.list.seas <- list(B = B, U = U, Q = Q, A = A, R = R, Z = Z, 
    C = C, c = c, x0 = x0, tinitx = 0)

Now we can fit the model:

m <- apply(dat.yr, 1, mean, na.rm = TRUE)
fit.seas <- MARSS(dat.yr, model = mod.list.seas, control = list(maxit = 500), 
    inits = list(A = matrix(m, ns, 1)))

Figure () shows the seasonal estimate plus mean for each station. This is \(z_i x_i + a_i\).

The estimated SWE at each station is \(E(y_i|y_{1:T})\). This is the estimate conditioned on all the data and includes the seasonal component plus the information from the data from other stations. Because we estimated a \(\mathbf{R}\) matrix with covariance, stations with data at time \(t\) help inform the value of stations without data at time \(t\). Only years up to 1990 are shown, but the model is fit to all years.

# this is the estimate of y conditioned on all the data
dd <- MARSShatyt(fit.seas)$ytT
rownames(dd) <- rownames(dat.yr)
colnames(dd) <- colnames(dat.yr)
ddd <- reshape2::melt(dd)
ddd$Var2 <- factor(ddd$Var2, levels = paste0(rep(1981:2013, each = 12), 
    "_", month.abb))
colnames(ddd) <- c("Station", "Year_Month", "SWE")
ddd <- ddd[order(ddd$Station, ddd$Year_Month), ]
ddd$Date <- swe.yr$Date
ddd$Year <- swe.yr$Year
ddd <- subset(ddd, Year < 1990)
p <- ggplot(data = ddd) + geom_line(aes(x = Date, y = SWE))
p <- p + geom_point(data = subset(swe.yr, Year < 1990), mapping = aes(x = Date, 
    y = SWE))
p + facet_wrap(~Station) + xlab("") + ylab("SWE")

Warning: Removed 1250 rows containing missing values (geom_point).